﻿namespace SpikingAnalyticsLib.ResponseFingerPrint

open System
open System.IO
open SpikingNeuronLib
open SpikingAnalyticsLib
open SpikingAnalyticsLib.GraphUtilities

//=================================================================================================

/// The set of windows that define a unique fingerprint for a given network and pattern
// Each neuron may have multiple windows
// Each window is defined by:
//  1. a list of temporal offsets within the frame
//  2. a list of spike counts for each offset
//  3. the total number of spikes in the frame
// Static methods exist to select the best window for each neuron or extract different data for each window
// e.g. windowSpikeCountRatio, frameSpikeFrequency, means and variance for each peak
// Notes:
// 1. The constructor returns RAW windows that should be filtered (typically with a final consistency threshold of >= 0.75)
// Use CreateFilteredWindowMap() to get the final window numbers for a specified consistency threshold
// Returning raw windows rather than filtered ones is more useful experimentally
//
// 2. The current algorithm returns MULTIPLE windows per neurons whereas the original algorithm
// returned the single best window (using a final consistency threshold of 0.75)
type WindowMap (frameProfile:FrameProfile, neuronWindowMap:Map<int, (int list * int list * int) list>) =

    // Extract the window raw data for the specified peaks found in the specified response histogram
    // part of NEW SelectWindows method
    static let GetRawWindowData (responseHistogram:Map<int, int>) (peaks:int list list) =
        peaks
        |> Seq.map (fun offsets ->
                let spikeCounts =
                    offsets
                    |> Seq.map (fun offset -> responseHistogram.[offset])
                    |> Seq.toList

                let totalSpikes =
                    responseHistogram
                    |> Map.fold (fun total offset count -> total + count) 0

                offsets, spikeCounts, totalSpikes
            )

    // Compute windowMinTime, windowMaxTime, window spike count ratio, frame spike frequency and peakiness
    static let ComputeSpikeRatios numFrames (rawData:int list * int list * int) =
        let offsets, spikeCounts, totalSpikes = rawData
        let windowMinTime, windowMaxTime = offsets.Head, offsets |> Seq.last

        let windowSpikes =
            spikeCounts
            |> Seq.sum

        // the spike count ratio: the proportion of profile spikes that occur in the selected window
        // (the conditional probability of window activation given the stimulus)
        // specificity
        let windowSpikeCountRatio = float windowSpikes / float totalSpikes

        // the stimulus response ratio: the ratio of the number of spikes in each window to the number of frames in the data
        // consistency
        let frameSpikeFrequency = float windowSpikes / float numFrames

        let peakiness = float windowSpikes / float (windowMaxTime - windowMinTime + 1)

        windowMinTime, windowMaxTime, windowSpikeCountRatio, frameSpikeFrequency, peakiness

    // Find consistent peaks using a small fixed-size window and remove overlaps between the selected offset ranges for each peak
    // Any windows found can be expanded using ExpandAndGroupWindows
    // part of NEW SelectWindows method
    static let FindActiveOffsets baseThreshold windowSize frameSize minSpikeCount (responseHistogram:Map<int, int>) =
        responseHistogram
        |> Map.toSeq
        |> Seq.windowed windowSize
        // count the total number of spikes
        |> Seq.map (fun window ->   // an array of (time, count)
                // as the window slides, get the spike count for each position
                let windowCount =
                    window
                    |> Seq.sumBy (fun (time, count) -> count)
                window, windowCount
            )
        // select windows that have high consistency across frames
        |> Seq.filter (fun (_, windowCount) -> windowCount >= minSpikeCount)
        // write the selected offset ranges to an array to remove overlaps
        |> Seq.fold (fun (active:bool []) (window, windowCount) ->
                // exclude offsets that have spike counts below the base threshold
                // note that this can leave isolated offsets with low overall counts
                // (therefore needs filtering with e.g. FilterByFrameSpikeFrequency)
                window
                |> Seq.iter (fun (offset, count) -> if count > baseThreshold then active.[offset] <- true)
                active
            ) (Array.zeroCreate frameSize)

    // Expand active offsets around each peak and group any overlapping offsets that belong to the same peak
    // part of NEW SelectWindows method
    static let ExpandAndGroupWindows frameSize baseThreshold (activeOffsets:bool []) (responseHistogram:Map<int, int>) =
        // expand trailing edges following each peak
        for i in 0..frameSize-2 do
            if activeOffsets.[i] then
                if activeOffsets.[i] <> activeOffsets.[i+1] then    // a trailing edge
                    let c1, c2 = responseHistogram.[i], responseHistogram.[i+1]
                    if c2 > baseThreshold && c2 <= c1 then
                        activeOffsets.[i+1] <- true

        // expand leading edges before each peak
        for i in frameSize-1..-1..1 do
            if activeOffsets.[i] then
                if activeOffsets.[i-1] <> activeOffsets.[i] then    // a leading edge
                    let c1, c2 = responseHistogram.[i-1], responseHistogram.[i]
                    if c1 > baseThreshold && c1 <= c2 then
                        activeOffsets.[i-1] <- true

        // group the expanded offset ranges
        activeOffsets
        |> Seq.mapi (fun index active -> index, active)
        |> Seq.fold (fun (groups:int list option list) (index, active) ->
                let currentWindow = groups.Head
                if active then
                    if currentWindow.IsNone then
                        Some([index]) :: groups.Tail
                    else
                        Some(index :: currentWindow.Value) :: groups.Tail
                else
                    if currentWindow.IsSome then
                        let orderedOffsets = List.rev currentWindow.Value
                        None :: Some(orderedOffsets) :: groups.Tail
                    else
                        groups
            ) [ None ]
        |> Seq.choose id
        |> Seq.toList

    // Select one or more windows for each neural response histogram (a window for each peak)
    // part of NEW SelectWindows method
    static let SelectWindowsImpl windowSize frameSize numFrames initialConsistencyThreshold baseThreshold excludeSet allProfiles =
        // compute the minimum spike count needed for initial fixed-sized windows to meet the required initial consistency estimate (window spikes per frame)
        // the initial consistency is typically kept low for the first thresholding step
        let minSpikeCount = int (initialConsistencyThreshold * float numFrames)

        allProfiles
        |> Map.toSeq
        |> Seq.filter (fun (nindex, _) -> not (Set.contains nindex excludeSet))
        |> Seq.map (fun (nindex, responseHistogram) ->
                // find peaks by finding temporal regions (windows) with high cumulative spike counts that exceed a threshold for consistency
                // get a list of offsets for each peak
                let peakOffsets =
                    // select consistent peaks using a fixed-size window
                    // remove overlaps between the offset ranges for each peak
                    let activeOffsets = FindActiveOffsets baseThreshold windowSize frameSize minSpikeCount responseHistogram
                    ExpandAndGroupWindows frameSize baseThreshold activeOffsets responseHistogram

                if peakOffsets.Length > 0 then
                    Some(nindex, GetRawWindowData responseHistogram peakOffsets)
                else
                    None
            )
        |> Seq.choose id

    // Get the window data for the specified profiles using pre-filtered windows from another frame profile
    // the pre-filtered windows define the allowed offsets for each response histogram (the regions of interest)
    // part of NEW SelectWindows method
    static let SelectWindowsPreFilteredImpl windowSize frameSize numFrames initialConsistencyThreshold baseThreshold (filterWindowMap:Map<int, #seq<int list * int list * int>>) allProfiles verbose =

        // Check that all pre-selected windows have at least partial overlap with the active offsets 
        let CheckAllPreSelectedWindowsHaveOverlap preSelectedWindowData (overlappingActiveOffsets:bool []) =
            preSelectedWindowData
            |> Seq.map (fun (allowedOffsetWindow, _, _) ->
                    // for each allowed window return true if there is at least some overlap with the active offsets
                    allowedOffsetWindow
                    |> Seq.exists (fun offset -> overlappingActiveOffsets.[offset])
                )
            // return true if ALL windows have at least some overlap
            |> Seq.forall (fun peakOverlapsAllowedOffsets -> peakOverlapsAllowedOffsets)

        // compute the minimum spike count needed for initial fixed-sized windows to meet the required initial consistency estimate (window spikes per frame)
        let minSpikeCount = int (initialConsistencyThreshold * float numFrames)

        let includeSet =
            filterWindowMap
            |> Map.toSeq
            |> Seq.map (fun (nindex, _) -> nindex)
            |> Set.ofSeq

        if verbose then printfn "No overlap with one or more windows for the following neurons:"

        allProfiles
        |> Map.toSeq
        |> Seq.filter (fun (nindex, _) -> Set.contains nindex includeSet)
        |> Seq.map (fun (nindex, responseHistogram) ->
                let peakOffsets =
                    // select consistent peaks using a fixed-size window
                    let activeOffsets =
                        // get the allowed offsets for the current response histogram
                        let allowedOffsets =
                            filterWindowMap.[nindex]
                            |> Seq.map (fun (offsets, _, _) -> offsets)
                            |> Seq.fold (fun (active:bool []) offsets ->
                                    offsets
                                    |> Seq.iter (fun offset -> active.[offset] <- true)
                                    active
                                ) (Array.zeroCreate frameSize)

                        // set active offsets that overlap with the allowed offsets
                        let overlappingOffsets =
                            FindActiveOffsets baseThreshold windowSize frameSize minSpikeCount responseHistogram
                            // filter the active offsets (only active offsets within the pre-filtered windows are allowed)
                            |> Seq.zip allowedOffsets
                            |> Seq.map (fun (insideWindow, active) -> insideWindow && active)
                            |> Seq.toArray

                        // warn if any pre-selected windows fail to overlap the active offsets
                        if not (CheckAllPreSelectedWindowsHaveOverlap filterWindowMap.[nindex] overlappingOffsets) then
                            if verbose then printfn "%d" nindex

                        overlappingOffsets

                    // expand and group the window for each peak
                    // note that this may extend outside of the allowed offsets
                    ExpandAndGroupWindows frameSize baseThreshold activeOffsets responseHistogram

                if peakOffsets.Length > 0 then
                    Some(nindex, GetRawWindowData responseHistogram peakOffsets)
                else
                    None
            )
        |> Seq.choose id

    // Create a new WindowMap
    // part of NEW SelectWindows method
    static let CreateWindowMap windowSize initialConsistencyThreshold baseThreshold (profile:FrameProfile) includeInhibitoryNeurons verbose =

        //let windowSize = profile.WindowSize
        //let numFrames = profile.NumberOfFrames
        //let responseHistograms = profile.ResponseHistograms

        let excludeSet =
            let patternSet =
                let pattern = profile.Pattern
                if pattern.IsSome then
                    pattern.Value.PatternNeurons
                    |> Set.ofSeq
                else
                    Set.empty
            let inhibitoryNeuronSet =
                if includeInhibitoryNeurons then
                    Set.empty
                else
                    Set.ofSeq (seq { for i in 800..999 do yield i })
            Set.union patternSet inhibitoryNeuronSet

        let neuronWindowMap =
            SelectWindowsImpl windowSize profile.FrameSize profile.NumberOfFrames initialConsistencyThreshold baseThreshold excludeSet profile.ResponseHistograms
            |> Seq.map (fun (nindex, windows) -> nindex, Seq.toList windows)
            |> Map.ofSeq

        if verbose then
            let windowCount =
                neuronWindowMap
                |> Map.fold (fun count _ windows -> count + windows.Length) 0
            printfn "Found %d windows on %d neurons" windowCount neuronWindowMap.Count

        neuronWindowMap

    // Take the midpoint of each window to create a sequence of firing events
    static let TransformWindowsToFiringEvents neuronWindowMap (frameProfile:FrameProfile) =

// ORIGINAL WindowMap version
//        neuronWindowMap
//        |> Map.toSeq
//        |> Seq.map (fun (neuronIndex, (min, max, linkProbability)) -> (min + max) / 2, neuronIndex, linkProbability)
//        |> Seq.toList

        neuronWindowMap
        |> Map.toSeq
        |> WindowMap.ExtractRatioData frameProfile.NumberOfFrames
        |> Seq.map (fun (neuronIndex, window) ->
                window
                |> Seq.map (fun (windowMinTime, windowMaxTime, windowSpikeCountRatio, frameSpikeFrequency, peakiness) ->
                        (windowMinTime + windowMaxTime) / 2, neuronIndex, windowSpikeCountRatio
                    )
            )
        |> Seq.concat
        |> Seq.toList

    static let GetFirstLayerWindowsFromNetwork (network:CrossbarNetwork) frameProfile synapticThreshold neuronWindowMap =
        let firstLevelEvents, _ =
            // take the midpoint of each window and create a sequence of firing events
            let windowFiringEvents = TransformWindowsToFiringEvents neuronWindowMap frameProfile
            network.GetLayeredFiringEvents(frameProfile.Pattern.Value, windowFiringEvents, synapticThreshold)

        // filter windows by the first-level firing event
        firstLevelEvents
        |> Seq.map (fun (time, nindex, _) ->
                let matchingWindowData =
                    neuronWindowMap.[nindex]
                    |> Seq.filter (fun (offsets, _, _) -> (Set.ofList offsets).Contains time)
                    |> Seq.toList
                if matchingWindowData.Length = 0 then
                    None
                else
                    if matchingWindowData.Length > 1 then
                        printfn "WARNING: event matched multiple windows for neuron %d" nindex
                    Some(nindex, matchingWindowData)
            )
        |> Seq.choose id
        |> Map.ofSeq

//---------------------------------------------------------------------------------------------------------

    // Get the overlap between temporal windows for this pair of profiles
    static let GetWindowSetIntersection (profile1:FrameProfile) (profile2:FrameProfile) temporalDifference =
        // get the windows for each profile
        let patternWindows =
            [ profile1; profile2 ]
            |> Seq.map (fun fp ->
                    let windowMap = new WindowMap(fp, false, false)
                    fp.PatternName, windowMap.AllWindows
                )
            |> Map.ofSeq

        let patternOverlap =
            // get the set of neurons that overlap for the patterns
            let overlappingNeurons =
                let patternWindowNeurons =
                    patternWindows
                    |> Map.toSeq
                    |> Seq.map (fun (patternName, windowMap) ->
                            windowMap
                            |> Map.toSeq
                            |> Seq.map (fun (neuronIndex, _) -> neuronIndex)
                            |> Set.ofSeq
                        )
                Set.intersectMany patternWindowNeurons

            // filter by overlapping neurons and invert key (key = neuronIndex; value = list of patternName, midTime)
            let neuronPatternMap =
                patternWindows
                |> Map.toSeq
                |> Seq.map (fun (patternName, windowMap) ->
                        windowMap
                        |> Map.toSeq
                        |> Seq.filter (fun (neuronIndex, _) -> overlappingNeurons.Contains neuronIndex)
                        |> Seq.map (fun (neuronIndex, windows) ->
                                windows
                                |> Seq.map (fun (offsets:int list, _, _) ->
                                        let windowMinTime, windowMaxTime = offsets.Head, offsets |> Seq.last
                                        neuronIndex, patternName, (windowMinTime + windowMaxTime) / 2
                                    )
                            )
                        |> Seq.concat
                    )
                |> Seq.concat
                |> Seq.fold (fun (state:Map<int, _>) (neuronIndex, patternName, midTime) ->
                        let data =
                            if state.ContainsKey neuronIndex then
                                state.[neuronIndex]
                            else
                                []
                        Map.add neuronIndex ((patternName, midTime) :: data) state
                    ) Map.empty

            // filter the data for each overlapping neuron
            // where the data consists of midTime values for each pattern
            // keep midTime values that are within 'temporalDifference' msec of each other
            // produces a sequence of patternName, point pairs where each point is (time, nindex)
            let filteredPoints =
                neuronPatternMap
                |> Map.toSeq
                |> Seq.map (fun (neuronIndex, matchedData) ->
                        // produce an array of lists, one for each pattern
                        // normally each list would have a single member (but could be more)
                        let firingTimeLists =
                            matchedData
                            |> Seq.fold (fun (state:Map<string, _>) (patternName, midTime) ->
                                    let data =
                                        if state.ContainsKey patternName then
                                            printfn "Warning: multiple values on %d %s" neuronIndex patternName
                                            state.[patternName]
                                        else
                                            []
                                    Map.add patternName (midTime :: data) state
                                ) Map.empty
                            |> Map.toSeq
                            |> Seq.toArray
                        // Warning: the following assumes that only two pattern profiles are present
                        // produce tuples from all pair-wise combinations from the two lists
                        seq {
                            let patternName1, midTimeList1 = firingTimeLists.[0]
                            let patternName2, midTimeList2 = firingTimeLists.[1]
                            for i in midTimeList1 do
                                for j in midTimeList2 do
                                    // pass through (patternName, (midTime, neuronIndex)) for midTime values with diff < 8
                                    if Math.Abs(i - j) <= temporalDifference then
                                        yield (patternName1, (i, neuronIndex))  // time, nindex
                                        yield (patternName2, (j, neuronIndex))
                        }
                    )
                |> Seq.concat

            // accumulate points for each pattern -> Map<patternName, point list>
            filteredPoints
            |> Seq.fold (fun (state:Map<string, _>) (patternName, point) ->
                    let pointList =
                        if state.ContainsKey patternName then
                            state.[patternName]
                        else
                            []
                    Map.add patternName (point :: pointList) state
                ) Map.empty

        patternOverlap

    // Find the optimum time range (window) of size windowSize for each neuron
    // Note that this currently doesn't allow for more than one window per neuron
    // ORIGINAL method
    static let SelectWindowsOriginal windowSize numFrames excludeSet allProfiles =
        allProfiles
        |> Map.toSeq
        |> Seq.filter (fun (nindex, _) -> not (Set.contains nindex excludeSet))
        |> Seq.map (fun (nindex, profile) ->
                // choose the window with the biggest window spike count
                let windowMinTime, windowMaxTime, windowSpikes =
                    profile
                    |> Map.toSeq
                    |> Seq.windowed windowSize
                    |> Seq.map (fun window ->   // an array of (time, count)
                            // get the start time, stop time and total spike count for the window
                            let minTime, _ = window.[0]
                            let maxTime, _ = window.[windowSize-1]
                            let windowCount =
                                window
                                |> Seq.sumBy (fun (time, count) -> count)
                            minTime, maxTime, windowCount
                        )
                    |> Seq.maxBy (fun (minTime, maxTime, windowCount) -> windowCount)

                let totalSpikes =
                    profile
                    |> Map.fold (fun total time count -> total + count) 0

                // the spike count ratio: the proportion of profile spikes that occur in the selected window
                // specificity
                let windowSpikeCountRatio = float windowSpikes / float totalSpikes

                // the stimulus response ratio: the ratio of the number of spikes in each window to the number of frames in the data
                // consistency
                let frameSpikeFrequency = float windowSpikes / float numFrames

                nindex, (windowMinTime, windowMaxTime, windowSpikeCountRatio, frameSpikeFrequency)
            )

    // Filter the selected neuron windows keeping windows in which the proportion of spikes in the window exceeds some threshold
    // ORIGINAL method
    static let FilterWindows initialThreshold frameFrequencyThreshold windows =
        windows
        |> Seq.filter (fun (_, (_, _, windowSpikeCountRatio, frameSpikeFrequency)) ->
                (windowSpikeCountRatio >= initialThreshold) && (frameSpikeFrequency >= frameFrequencyThreshold)
            )
        |> Seq.map (fun (nindex, (minTime, maxTime, windowSpikeCountRatio, frameSpikeFrequency)) ->
                nindex, (minTime, maxTime, windowSpikeCountRatio)
            )
        |> Map.ofSeq

    static let CreateWindowMapOriginal frameFrequencyThreshold (profile:FrameProfile) includeInhibitoryNeurons verbose =

        let windowSize = profile.WindowSize
        let numFrames = profile.NumberOfFrames
        let responseHistograms = profile.ResponseHistograms

        let excludeSet =
            let patternSet =
                let pattern = profile.Pattern
                if pattern.IsSome then
                    pattern.Value.PatternNeurons
                    |> Set.ofSeq
                else
                    Set.empty
            let inhibitoryNeuronSet =
                if includeInhibitoryNeurons then
                    Set.empty
                else
                    Set.ofSeq (seq { for i in 800..999 do yield i })
            Set.union patternSet inhibitoryNeuronSet

        let neuronWindowMap =
            // select windows on the selected neurons
            let windows =
                SelectWindowsOriginal windowSize numFrames excludeSet responseHistograms
                |> Seq.toList

            // create an initial threshold (to reduce the amount of processing)
            // set the initial threshold to the average
            let initialThreshold =
                windows
                |> Seq.map (fun (nindex, (_, _, spikeCountRatio, _)) -> spikeCountRatio)
                |> Seq.average

            if verbose then
                printfn "Average spike count ratio: %f" initialThreshold

            // filter the windows
            windows
            |> FilterWindows initialThreshold frameFrequencyThreshold

            // CHANGES for compatibility with NEW SelectWindows method
            |> Map.map (fun nindex (minTime, maxTime, _) ->
                    let peakOffsets = [ for index in minTime..maxTime -> index ]
                    let responseHistogram = responseHistograms.[nindex]

                    let spikeCounts =
                        peakOffsets
                        |> Seq.map (fun offset -> responseHistogram.[offset])
                        |> Seq.toList

                    let totalSpikes =
                        responseHistogram
                        |> Map.fold (fun total offset count -> total + count) 0

                    [ peakOffsets, spikeCounts, totalSpikes ]
                )

        neuronWindowMap

//---------------------------------------------------------------------------------------------------------
// Processing methods

    // Extract the window ratio data for the specified raw data generated by GetWindowRawData
    static let ExtractRatioDataImpl numFrames (rawData:#seq<(int * #seq<int list * int list * int>)>) =
            rawData
            |> Seq.map (fun (nindex, windows) ->
                let ratioData =
                    windows
                    |> Seq.map (fun windowData ->
                            ComputeSpikeRatios numFrames windowData
                        )

                nindex, ratioData
            )

    // Extract the window mean and variance data for the specified raw data generated by GetWindowRawData
    static let ExtractMeanVarianceDataImpl (rawData:#seq<int * #seq<int list * int list * int>>) =
        rawData
        |> Seq.map (fun (nindex, windowData) ->
            let meanVarianceData =
                windowData
                |> Seq.map (fun (offsets, spikeCounts, totalSpikes) ->
                        let flattenedOffsets =
                            if offsets.Length > 1 then
                                Seq.zip offsets spikeCounts
                                |> Seq.map (fun (offset, spikeCount) ->
                                        // repeat each offset spikeCount times
                                        Seq.init spikeCount (fun index -> float offset)
                                    )
                                |> Seq.concat
                                |> Seq.toList
                            else
                                offsets
                                |> List.map (fun offset -> float offset)

                        let mean =
                            flattenedOffsets
                            |> List.average

                        let variance =
                            let sumSquareDifferences =
                                flattenedOffsets
                                |> Seq.sumBy (fun offset ->
                                        let diff = offset - mean
                                        diff * diff
                                    )
                            // population variance
                            let N = float flattenedOffsets.Length
                            sumSquareDifferences / N

                        let windowSpikes =
                            spikeCounts
                            |> Seq.sum

                        mean, variance, windowSpikes, totalSpikes
                    )

            nindex, meanVarianceData
        )

    static let ChooseBestRawWindowsImpl (numFrames:int) (rawData:#seq<int * (int list * int list * int) list>) =
        rawData
        |> Seq.map (fun (nindex, windowData) ->
            let offsets, spikeCounts, totalSpikes, frameSpikeFrequency =
                windowData
                |> Seq.map (fun (offsets, spikeCounts, totalSpikes) ->

                        let windowSpikes =
                            spikeCounts
                            |> Seq.sum

                        let frameSpikeFrequency = float windowSpikes / float numFrames

                        offsets, spikeCounts, totalSpikes, frameSpikeFrequency
                    )
                |> Seq.maxBy (fun (offsets, spikeCounts, totalSpikes, frameSpikeFrequency) -> frameSpikeFrequency)

            nindex, (offsets, spikeCounts, totalSpikes)
        )

    // Choose a single window for each neural response histogram (the window that has the best windowSpikeCountRatio)
    static let ChooseBestRatioWindowsImpl (ratioData:#seq<int * #seq<int * int * float * float * float>>) =
        ratioData
        |> Seq.map (fun (nindex, windows) ->
                let bestWindow =
                    windows
                    |> Seq.maxBy (fun (windowMinTime, windowMaxTime, windowSpikeCountRatio, frameSpikeFrequency, peakiness) -> windowSpikeCountRatio)
                nindex, bestWindow
            )

    // Choose a single window for each neural response histogram (the window that has the smallest variance)
    static let ChooseBestVarianceWindowsImpl (meanVarianceData:seq<int * seq<float * float * int * int>>) =
        meanVarianceData
        |> Seq.map (fun (nindex, windows) ->
                let bestWindow =
                    windows
                    |> Seq.minBy (fun (mean, variance, windowSpikes, totalSpikes) -> variance)
                nindex, bestWindow
            )

    static let CreateMeanVarianceMap windows =
        windows
        |> ExtractMeanVarianceDataImpl
        |> ChooseBestVarianceWindowsImpl
        |> Map.ofSeq

//---------------------------------------------------------------------------------------------------------
// Constructors

    /// Select windows using the specified frame profile
    new (frameProfile, includeInhibitoryNeurons) =

        let initialConsistencyThreshold = 0.5
        let baseThreshold = 4
        let windowSize = 8

        let neuronWindowMap = CreateWindowMap windowSize initialConsistencyThreshold baseThreshold frameProfile includeInhibitoryNeurons false
        new WindowMap(frameProfile, neuronWindowMap)

    /// Select windows using the specified frame profile
    new (frameProfile, includeInhibitoryNeurons, verbose) =

        let initialConsistencyThreshold = 0.5
        let baseThreshold = 4
        let windowSize = 8

        let neuronWindowMap = CreateWindowMap windowSize initialConsistencyThreshold baseThreshold frameProfile includeInhibitoryNeurons verbose
        new WindowMap(frameProfile, neuronWindowMap)

    /// Select windows using the specified frame profile
    new (initialConsistencyThreshold, baseThreshold, windowSize, frameProfile:FrameProfile, includeInhibitoryNeurons, verbose) =

        assert (initialConsistencyThreshold > 0.0 && initialConsistencyThreshold <= 1.0)
        assert (baseThreshold >= 0)
        assert (windowSize > 0 && windowSize < frameProfile.FrameSize)

        let neuronWindowMap = CreateWindowMap windowSize initialConsistencyThreshold baseThreshold frameProfile includeInhibitoryNeurons verbose
        new WindowMap(frameProfile, neuronWindowMap)

//---------------------------------------------------------------------------------------------------------
// Properties

    /// The name of the pattern that matches these windows
    member this.PatternName = frameProfile.PatternName

    /// The frame profile for this window map and pattern
    member this.FrameProfile = frameProfile

    /// All selected windows keyed by neuron
    // there may be multiple windows per neurons
    member this.AllWindows = neuronWindowMap

    /// All selected windows as a list keyed by neuron
    // there may be multiple windows per neurons
    member this.AllWindowsList =
        neuronWindowMap
        |> Map.toList

    /// All selected windows as a list keyed by neuron
    // there may be multiple windows per neurons
    member this.AllWindowsSequence =
        neuronWindowMap
        |> Map.toSeq

    /// All selected windows keyed first by neuron and then window index
    // The window index is assigned by the rank order of each window peak in the list of windows for each neuron
    member this.AllWindowsIndexed =
        neuronWindowMap
        |> Map.map (fun nindex windows ->
                windows
                |> Seq.mapi (fun index window -> index, window)
                |> Map.ofSeq
            )

    /// The single best window for each neuron (the window that produces the best windowSpikeCountRatio)
    member this.BestRatioWindows =
        neuronWindowMap
        |> Map.toSeq
        |> WindowMap.ExtractRatioData frameProfile.NumberOfFrames
        |> WindowMap.ChooseBestRatioWindows
        |> Map.ofSeq

    /// The single best window for each neuron (the window that produces the smallest variance)
    member this.SmallestVarianceWindows =
        neuronWindowMap
        |> Map.toSeq
        |> WindowMap.ExtractMeanVarianceData
        |> WindowMap.ChooseBestVarianceWindows
        |> Map.ofSeq

    /// Get the neurons for the selected windows for this pattern
    member this.WindowNeurons
        with get() =
            neuronWindowMap
            |> Map.toSeq
            |> Seq.map (fun (nindex, _) -> nindex)
            |> Seq.toList

    /// Write a brief summary to the console
    member this.Summary =
        let peakCount =
            neuronWindowMap
            |> Map.toSeq
            |> Seq.map (fun (nindex, windows:_ list) ->
                    windows.Length
                )
            |> Seq.max
        printfn "Neurons with windows: %d\nMaximum windows per neuron: %d" (neuronWindowMap.Count) peakCount

    /// Get a function to extract mean and variance data from the raw data returned by GetWindowData
    static member ExtractMeanVarianceData =
        ExtractMeanVarianceDataImpl

    /// Get a function to extract ratio data from the raw data returned by GetWindowData
    static member ExtractRatioData =
        ExtractRatioDataImpl

    /// Get a function to choose the best windows from the data returned by ExtractMeanVarianceData
    static member ChooseBestVarianceWindows =
        ChooseBestVarianceWindowsImpl

    /// Get a function to choose the best windows from the data returned by ExtractRatioData
    static member ChooseBestRatioWindows =
        ChooseBestRatioWindowsImpl

    /// Get a function to choose the best windows from the raw data (for all neurons)
    static member ChooseBestRawWindows =
        ChooseBestRawWindowsImpl

    /// Get a function to choose the best window from the raw data for a single neuron
    static member ChooseBestRawWindow =
        let ChooseBestRawWindow (numFrames:int) (windowData:(int list * int list * int) list) =
            let offsets, spikeCounts, totalSpikes, frameSpikeFrequency =
                windowData
                |> Seq.map (fun (offsets, spikeCounts, totalSpikes) ->

                        let windowSpikes =
                            spikeCounts
                            |> Seq.sum

                        let frameSpikeFrequency = float windowSpikes / float numFrames

                        offsets, spikeCounts, totalSpikes, frameSpikeFrequency
                    )
                |> Seq.maxBy (fun (offsets, spikeCounts, totalSpikes, frameSpikeFrequency) -> frameSpikeFrequency)

            offsets, spikeCounts, totalSpikes

        ChooseBestRawWindow

    /// Get a function to select windows that overlap the single window from the reference data, using the raw data for a single neuron
    // Any window that fails to overlap the reference window is rejected
    static member SelectOverlappingWindows =
        let SelectOverlappingWindows (referenceData:int list * int list * int) (windowData:(int list * int list * int) list) =
            let IsOverlap offsets1 offsets2 =

                let allOffsets =
                    [offsets1; offsets2]
                    |> Seq.concat
                    |> Seq.toList

                let min = Seq.min allOffsets
                let max = Seq.max allOffsets

                let active = Array.zeroCreate (max+1)

                // set active offsets
                allOffsets
                |> Seq.iter (fun offset -> active.[offset] <- true)

                // check that no inactive offsets between min and max
                active
                |> Seq.mapi (fun index active -> index, active)
                |> Seq.forall (fun (index, active) -> if index >= min then active else true)

            let referenceOffsets, _, _ = referenceData
            windowData
            |> Seq.filter (fun (offsets, _, _) -> IsOverlap referenceOffsets offsets)
            |> Seq.toList

        SelectOverlappingWindows

    /// Get a function to count the number of spikes in the central mass of an isolated spike count peak
    // Uses a sliding window to find the maximum spike count in the spike count list
    static member CountCentralMass =
        let CountCentralMass windowSize (spikeCounts:int list) =
            if spikeCounts.Length >= windowSize then
                spikeCounts
                |> Seq.windowed windowSize
                |> Seq.map (fun (window:int []) ->   // an array of (time, count)
                        Seq.sum window
                    )
                |> Seq.max
            else
                Seq.sum spikeCounts

        CountCentralMass

    /// Get a function to filter the raw data returned from GetWindowData by the window spikes per frame (consistency)
    static member FilterByFrameSpikeFrequency =
        let FilterData numFrames threshold (rawData:#seq<int * (int list * int list * int) list>) =
            rawData
            |> Seq.map (fun (nindex, windowData) ->
                    let filteredWindowData =
                        windowData
                        |> List.filter (fun (offsets, spikeCounts, totalSpikes) ->
                                let windowSpikes =
                                    spikeCounts
                                    |> Seq.sum

                                let frameSpikeFrequency = float windowSpikes / float numFrames   // consistency: window spikes per frame

                                frameSpikeFrequency > threshold
                            )
                    nindex, filteredWindowData
                )
            |> Seq.filter (fun (nindex, windowData) -> windowData.Length > 0)

        FilterData

    /// Get a function to filter the raw data returned from GetWindowData by the window spikes per total spikes (selectivity)
    static member FilterByWindowSpikeCountRatio =
        let FilterData threshold (rawData:#seq<int * (int list * int list * int) list>) =
            rawData
            |> Seq.map (fun (nindex, windowData) ->
                    let filteredWindowData =
                        windowData
                        |> List.filter (fun (offsets, spikeCounts, totalSpikes) ->
                                let windowSpikes =
                                    spikeCounts
                                    |> Seq.sum

                                let windowSpikeCountRatio = float windowSpikes / float totalSpikes

                                windowSpikeCountRatio > threshold
                            )
                    nindex, filteredWindowData
                )
            |> Seq.filter (fun (nindex, windowData) -> windowData.Length > 0)

        FilterData

//---------------------------------------------------------------------------------------------------------
// Methods

    /// Filter the firing events that produce window activation (assuming a single frame of data)
    /// The jitter value may be zero
    member this.GetWindowActivationEvents(firingEvents, jitter) =

        let AllowJitter (offsets:int list) jitter =
            offsets
            |> Seq.map (fun offset ->
                set [ for i in offset - jitter..offset + jitter -> i ]
            )
            |> Set.unionMany

        // construct a frame offset map using the specified firing events
        // assuming a single frame of data i.e. firing times are not folded into the frame
        let frameOffsetMap =
            firingEvents
            |> Seq.fold (fun (state:Map<int,int list>) (time, nindex, probability) ->
                    // produce a list of reponse times for each neuron in the current frame
                    let times =
                        if state.ContainsKey(nindex) then
                            state.[nindex]
                        else
                            []
                    Map.add nindex (time :: times) state
                ) Map.empty

        // get each firing event that produced a window activation
        neuronWindowMap
        |> Map.toSeq
        |> Seq.map (fun (neuronIndex, windows) ->
                windows
                |> Seq.map (fun (offsets, _, _) ->
                        let windowSpikeTimes =
                            if frameOffsetMap.ContainsKey neuronIndex then
                                let times = Set.ofList frameOffsetMap.[neuronIndex]
                                let window = AllowJitter offsets jitter
                                Set.intersect window times
                                |> Seq.toList
                            else
                                []
                        windowSpikeTimes, neuronIndex
                    )
            )
        |> Seq.concat
        |> Seq.filter (fun (windowSpikeTimes, _) -> not (List.isEmpty windowSpikeTimes))
        |> Seq.map (fun (windowSpikeTimes, neuronIndex) ->
                windowSpikeTimes
                |> Seq.map (fun time -> time, neuronIndex)
            )
        |> Seq.concat
        |> Seq.toList

    /// Transform windows into firing events by taking the mid-point of each window
    member this.GetWindowsAsFiringEvents() = TransformWindowsToFiringEvents neuronWindowMap frameProfile

    /// Transform windows into firing events and merge in the stimulation pattern events
    // Note that the pattern events are the times of current onset and not the actual firing times
    // of the stimulated neurons; therefore add an appropriate offset to account for firing latency
    // This applies to the pattern events only i.e. not to the windowMap windows
    member this.GetWindowsAsFiringEvents(pattern:Pattern, ?latency) =
        let firingDelay = if latency.IsSome then latency.Value else 0
        pattern.GetFiringEvents()
        |> Seq.map (fun ev -> ev.Time + firingDelay, ev.NeuronIndex, 1.0)
        |> Seq.toList
        |> List.append (TransformWindowsToFiringEvents neuronWindowMap frameProfile)

    /// Use the window centre-points to generate a description of the PNG structure
    member this.GetPNGDescriptor(network:CrossbarNetwork, pattern:Pattern) =
        let fingerprintFiringEvents = this.GetWindowsAsFiringEvents(pattern)
        network.GetPNGDescriptor(fingerprintFiringEvents, 0.0, network.MaximumAxonalDelay)

    /// Compute the longest path through the PNG structure generated using the window centre-points
    // The specified triggering pattern is included in the PNG structure
    member this.GetLongestPath(network:CrossbarNetwork, trigger:Pattern) =
        // get a descriptor for the network structure
        let descriptor = this.GetPNGDescriptor(network, trigger)
        descriptor.GetLongestPath(trigger)

    /// Check if a single presentation of a stimulus derived from the specified pattern can activate all windows
    /// Returns the number of activated windows and the total number of windows
    member this.CheckActivation(network:CrossbarNetwork, pattern:Pattern, jitter) =
        // run the network with no background to avoid priming effect (smoothes v and u)
        network.Run(10, None, 0, false)

        // present a single trigger stimulus and collect firing events
        let stimulus = Stimulus.Create(1, pattern) :> IStimulus
        network.Run(1, Some(stimulus), 0, false)
        let firingEvents = network.EventTuples
        let activatedWindows = this.GetWindowActivationEvents(firingEvents, jitter)
        activatedWindows.Length, neuronWindowMap.Count

    /// Get the subset of firing events that correpond to a fingerprint window
    /// Which of the specified firing event matches a fingerprint window?
    // MG similar to GetWindowActivationEvents(firingEvents, jitter)
    member this.GetMatchedFiringEvents(firingEvents:FiringEvent list, pattern:Pattern, jitter) =

        // create a firing time map of fingerprint windows
        // add back the pattern events (pattern events are not profiled and therefore don't generate windows)
        let fingerprintNeuronTimeMap =
            this.GetWindowsAsFiringEvents(pattern)
            |> Seq.fold (fun (state:Map<int, int list>) (time, nindex, _) ->
                let timeList =
                    if state.ContainsKey nindex then
                        state.[nindex]
                    else
                        []
                Map.add nindex (time :: timeList) state
            ) Map.empty

        firingEvents
        |> Seq.map (fun ev -> ev.Time, ev.NeuronIndex)
        |> Seq.filter (fun (time, nindex) ->
            if fingerprintNeuronTimeMap.ContainsKey nindex then
                fingerprintNeuronTimeMap.[nindex]
                |> Seq.exists (fun _time -> time >= _time - jitter && time <= _time + jitter)
            else
                false
        )
        |> Set.ofSeq

    /// Generate firing events for each window and group them according to the underlying network structure
    // Group firing events into first layer (directly connected to the input stimulus firing events)
    // and the remainder (not directly connected to stimulus firing events)
    // The profile data must have been produced from the specified network
    member this.GetLayeredFiringEvents(network:CrossbarNetwork, synapticThreshold) =

        // take the midpoint of each window and create a sequence of firing events
        let windowFiringEvents = this.GetWindowsAsFiringEvents()
        network.GetLayeredFiringEvents(this.FrameProfile.Pattern.Value, windowFiringEvents, synapticThreshold)

    /// Create a new WindowMap by filtering this WindowMap to keep only first layer windows
    // i.e. windows whose firing events are directly connected to the input stimulus firing events
    member this.GetFirstLayerWindows(network, synapticThreshold) =
        let filteredNeuronWindowMap = GetFirstLayerWindowsFromNetwork network frameProfile synapticThreshold neuronWindowMap
        new WindowMap(frameProfile, filteredNeuronWindowMap)

//---------------------------------------------------------------------------------------------------------
// Static Methods

    /// Generate the window spike ratios for the specified window data
    static member ComputeWindowSpikeRatios(numFrames, rawData:int list * int list * int) =
        ComputeSpikeRatios numFrames rawData

    /// For each neuron find the optimum time range (window) that maximises the proportion of spikes in the window
    // Provided for compatibility with older APIs
    static member SelectBestWindowsOriginal(windowSize, numFrames, initialThreshold, frameFrequencyThreshold, excludeSet, allProfiles) =
        SelectWindowsOriginal windowSize numFrames excludeSet allProfiles
        |> FilterWindows initialThreshold frameFrequencyThreshold

    /// Create a new WindowMap by selecting windows using the ORIGINAL method and the specified frame profile
    static member CreateWithOriginalMethod(frameProfile, includeInhibitoryNeurons, verbose) =

        // final threshold value
        let frameFrequencyThreshold = 0.75

        let neuronWindowMap = CreateWindowMapOriginal frameFrequencyThreshold frameProfile includeInhibitoryNeurons verbose
        new WindowMap(frameProfile, neuronWindowMap)

    /// Create a new WindowMap by selecting windows using the ORIGINAL method and the specified frame profile
    static member CreateWithOriginalMethod(frameFrequencyThreshold, frameProfile, includeInhibitoryNeurons, verbose) =

        assert (frameFrequencyThreshold > 0.0 && frameFrequencyThreshold <= 1.0)

        let neuronWindowMap = CreateWindowMapOriginal frameFrequencyThreshold frameProfile includeInhibitoryNeurons verbose
        new WindowMap(frameProfile, neuronWindowMap)

    /// Get the overlap between temporal windows for this pair of profiles
    // allow the centre of each window to differ by up to 'temporalDifference' msecs
    static member GetWindowOverlap(profile1:FrameProfile, profile2:FrameProfile, temporalDifference) =
        GetWindowSetIntersection profile1 profile2 temporalDifference

    /// Transform the windows for the specified profile into firing events by taking the mid-point of each window
    static member GetWindowsAsFiringEvents(profile:FrameProfile, includeInhibitoryNeurons:bool) =

        let windowMap = new WindowMap(profile, includeInhibitoryNeurons, false)

        // take the midpoint of each window and create a sequence of firing events
        windowMap.GetWindowsAsFiringEvents()

    /// Transform the windows for the specified profile into firing events
    /// and group them into first layer (directly connected to the input stimulus firing events)
    /// and the remainder (not directly connected to stimulus firing events)
    // the profile data must have been produced from the specified network
    static member LayerWindowsAsFiringEvents(profile, network:CrossbarNetwork, synapticThreshold, includeInhibitoryNeurons) =

        let windowMap = new WindowMap(profile, includeInhibitoryNeurons, false)
        windowMap.GetLayeredFiringEvents(network, synapticThreshold)

    /// Get the window data (offsets, spike counts, total spikes) for each peak in the specified response histogram
    static member GetWindowData(responseHistogram:Map<int, int>, peakOffsets:int list list) =
        GetRawWindowData responseHistogram peakOffsets
        |> Seq.toList

    /// Find all windows for the specified profile (the response histogram for a single neuron)
    // May return multiple windows if more than one window meets the threshold criteria
    // This method is mainly for algorithm testing purposes as it calls into the window selection code
    static member SelectAllProfileWindows(responseHistogram:Map<int, Map<int, int>>, frameSize, numberOfFrames, windowSize, initialConsistencyThreshold, baseThreshold) =
        SelectWindowsImpl windowSize frameSize numberOfFrames initialConsistencyThreshold baseThreshold (Set.empty) responseHistogram
        |> Seq.map (fun (nindex, windows) -> nindex, Seq.toList windows)
        |> Seq.toList

    /// Find peaks over the specified threshold for the specified response histogram
    /// Write the offsets and associated counts to the console
    // This method is mainly for algorithm testing purposes as it calls into the window selection code
    static member ShowPeakOffsets(profile:FrameProfile, neuronIndex, initialConsistencyThreshold, baseThreshold) =
        let peakOffsets =
            let windowSize, frameSize, numFrames = profile.WindowSize, profile.FrameSize, profile.NumberOfFrames
            let responseHistogram = profile.ResponseHistograms.[neuronIndex]
            let minSpikeCount = int (initialConsistencyThreshold * float numFrames)     // window spikes per frame
            let activeOffsets = FindActiveOffsets baseThreshold windowSize frameSize minSpikeCount responseHistogram
            ExpandAndGroupWindows frameSize baseThreshold activeOffsets responseHistogram

        WindowMap.ShowPeakOffsets(profile.ResponseHistograms.[neuronIndex], peakOffsets)

    /// Write the specified peaks and associated counts to the console
    static member ShowPeakOffsets(responseHistogram:Map<int, int>, peakOffsets:int list list) =
        printfn "Number of peaks: %d" peakOffsets.Length
        peakOffsets
        |> Seq.map (fun peak ->
                let counts =
                    peak
                    |> Seq.map (fun offset -> responseHistogram.[offset])
                    |> Seq.toList
                peak, counts
            )
        |> Seq.iteri (fun index (offsets, counts) ->
            printfn "\nPeak %d" (index + 1)
            printfn "Offsets: %s" (System.String.Join(" ", offsets))
            printfn "Counts: %s" (System.String.Join(" ", counts))
        )

    /// Create a WindowMap where the windows have been filtered according to the specified consistency threshold
    // A typical final consistency threshold is >= 0.75
    static member CreateFilteredWindowMap(profile:FrameProfile, finalConsistencyThreshold, includeInhibitory, verbose) =
        let filteredWindows =
            let windowMap = new WindowMap(profile, includeInhibitory, verbose)
            windowMap.AllWindowsSequence
            |> WindowMap.FilterByFrameSpikeFrequency profile.NumberOfFrames finalConsistencyThreshold
            |> Map.ofSeq

        new WindowMap(profile, filteredWindows)

    /// Create a new WindowMap for the specified profile where window selection is filtered using window data from an existing WindowMap that was created from another profile
    static member CreateFilteredWindowMap(profile:FrameProfile, windowSize, initialConsistencyThreshold, finalConsistencyThreshold, baseThreshold, windowMap:WindowMap, verbose) =
        let filteredWindowMap =
            SelectWindowsPreFilteredImpl windowSize profile.FrameSize profile.NumberOfFrames initialConsistencyThreshold baseThreshold (windowMap.AllWindows) profile.ResponseHistograms verbose
            |> Seq.map (fun (nindex, windows) -> nindex, Seq.toList windows)
            |> WindowMap.FilterByFrameSpikeFrequency (profile.NumberOfFrames) finalConsistencyThreshold
            |> Map.ofSeq
        new WindowMap(profile, filteredWindowMap)

    /// Find the corresponding windows from the reference window map in the comparison profile
    // The comparison WindowMap is filtered so that only overlapping windows are included
    static member CompareProfileWindows(initialConsistencyThreshold, finalConsistencyThreshold, baseThreshold, windowSize, referenceWindowMap:WindowMap, comparisonProfile:FrameProfile, verbose) =
        let comparisonWindowMap = WindowMap.CreateFilteredWindowMap(comparisonProfile, windowSize, initialConsistencyThreshold, finalConsistencyThreshold, baseThreshold, referenceWindowMap, verbose)

        referenceWindowMap, comparisonWindowMap

    /// Find the corresponding windows from the reference window map in the comparison profile
    // The comparison WindowMap is filtered so that only overlapping windows are included
    static member CompareProfileWindows(referenceWindowMap:WindowMap, comparisonProfile:FrameProfile, verbose) =
        let initialConsistencyThreshold = 0.5
        let finalConsistencyThreshold = 0.75
        let baseThreshold = 4                       // minimum spike counts to be included at the edge of a peak
        let windowSize = 4                          // size of initial fixed-sized window scan

        WindowMap.CompareProfileWindows(initialConsistencyThreshold, finalConsistencyThreshold, baseThreshold, windowSize, referenceWindowMap, comparisonProfile, verbose)

    /// Select only first-level windows from the reference profile and then find the corresponding windows in the comparison profile
    // First-level windows are those collecting firing events that are directly connected to stimulus firing events
    // The comparison WindowMap is filtered so that only overlapping windows are included
    static member CompareProfileWindows(initialConsistencyThreshold, finalConsistencyThreshold, baseThreshold, windowSize, synapticThreshold, network:CrossbarNetwork, referenceProfile:FrameProfile, comparisonProfile:FrameProfile, verbose) =

        let referenceWindowMap =
            let referenceWindows =
                CreateWindowMap windowSize initialConsistencyThreshold baseThreshold referenceProfile false true
                |> GetFirstLayerWindowsFromNetwork network referenceProfile synapticThreshold
                |> Map.toSeq
                |> WindowMap.FilterByFrameSpikeFrequency (referenceProfile.NumberOfFrames) finalConsistencyThreshold
                |> Map.ofSeq

            new WindowMap(referenceProfile, referenceWindows)

        WindowMap.CompareProfileWindows(initialConsistencyThreshold, finalConsistencyThreshold, baseThreshold, windowSize, referenceWindowMap, comparisonProfile, verbose)

    /// Select only first-level windows from the reference profile and then find the corresponding windows in the comparison profile
    // First-level windows are those collecting firing events that are directly connected to stimulus firing events
    // The comparison WindowMap is filtered so that only overlapping windows are included
    static member CompareProfileWindows(network:CrossbarNetwork, referenceProfile:FrameProfile, comparisonProfile:FrameProfile, verbose) =

        let initialConsistencyThreshold = 0.5
        let finalConsistencyThreshold = 0.75
        let baseThreshold = 4                       // minimum spike counts to be included at the edge of a peak
        let synapticThreshold = 9.9
        let windowSize = 4                          // size of initial fixed-sized window scan

        WindowMap.CompareProfileWindows(initialConsistencyThreshold, finalConsistencyThreshold, baseThreshold, windowSize, synapticThreshold, network, referenceProfile, comparisonProfile, verbose)

    /// Select windows from the reference profile and then find the corresponding windows in the comparison profile
    // The comparison WindowMap is filtered so that only overlapping windows are included
    static member CompareProfileWindows(initialConsistencyThreshold, finalConsistencyThreshold, baseThreshold, windowSize, referenceProfile:FrameProfile, comparisonProfile:FrameProfile, verbose) =

        let referenceWindowMap =
            let referenceWindows =
                CreateWindowMap windowSize initialConsistencyThreshold baseThreshold referenceProfile false verbose
                |> Map.toSeq
                |> WindowMap.FilterByFrameSpikeFrequency (referenceProfile.NumberOfFrames) finalConsistencyThreshold
                |> Map.ofSeq

            new WindowMap(referenceProfile, referenceWindows)

        WindowMap.CompareProfileWindows(initialConsistencyThreshold, finalConsistencyThreshold, baseThreshold, windowSize, referenceWindowMap, comparisonProfile, verbose)

    /// Select windows from the reference profile and then find the corresponding windows in the comparison profile
    // The comparison WindowMap is filtered so that only overlapping windows are included
    static member CompareProfileWindows(referenceProfile:FrameProfile, comparisonProfile:FrameProfile, verbose) =

        let initialConsistencyThreshold = 0.5
        let finalConsistencyThreshold = 0.75
        let baseThreshold = 4                       // minimum spike counts to be included at the edge of a peak
        let windowSize = 4                          // size of initial fixed-sized window scan

        WindowMap.CompareProfileWindows(initialConsistencyThreshold, finalConsistencyThreshold, baseThreshold, windowSize, referenceProfile, comparisonProfile, verbose)

    /// Get the profile comparison data generated by CompareProfileWindows as a mean-variance comparison
    static member GetMeanVarianceData(referenceWindowMap:WindowMap, comparisonWindowMap:WindowMap) =
        let referenceMap, comparisonMap = CreateMeanVarianceMap referenceWindowMap.AllWindowsSequence, CreateMeanVarianceMap comparisonWindowMap.AllWindowsSequence

        // merge the two data sets
        // Note that Deedle Frame.ofRecords seems to be limited to 8 columns (additional columns more than that become a single tuple column)
        referenceMap
        |> Map.toSeq
        |> Seq.map (fun (nindex, referenceData) ->
                if comparisonMap.ContainsKey nindex then
                    Some(nindex, referenceData, comparisonMap.[nindex])
                else
                    None
            )
        |> Seq.choose id
        |> Seq.toList

    /// Save the profile comparison data generated by CompareProfileWindows as a mean-variance comparison
    static member SaveMeanVarianceData(filePath, referenceWindowMap:WindowMap, comparisonWindowMap:WindowMap) =

        let referenceMap, comparisonMap = CreateMeanVarianceMap referenceWindowMap.AllWindowsSequence, CreateMeanVarianceMap comparisonWindowMap.AllWindowsSequence

        // merge the two data sets
        let result =
            referenceMap
            |> Map.toSeq
            |> Seq.map (fun (nindex, referenceData) ->
                    if comparisonMap.ContainsKey nindex then
                        Some(nindex, referenceData, comparisonMap.[nindex])
                    else
                        None
                )
            |> Seq.choose id
            |> Seq.map (fun (nindex, (mean1, variance1, windowSpikes1, totalSpikes1), (mean2, variance2, windowSpikes2, totalSpikes2)) ->
                    seq {
                        // NO meta, then WITH meta
                        yield sprintf "%d %d %f %f %d %d" nindex 0 mean1 variance1 windowSpikes1 totalSpikes1
                        yield sprintf "%d %d %f %f %d %d" nindex 1 mean2 variance2 windowSpikes2 totalSpikes2
                    }
                )
            |> Seq.concat
            |> Seq.toList

        let header = "nindex withMeta mean variance windowSpikes totalSpikes"
        File.WriteAllLines(filePath, header :: result)

    /// Save the profile comparison data generated by CompareProfileWindows as a comparison of raw peak value data
    static member SaveRawPeakValueData(filePath, numFrames, referenceWindowMap:WindowMap, comparisonWindowMap:WindowMap) =

        let referenceBestWindows =
            referenceWindowMap.AllWindowsSequence
            |> WindowMap.ChooseBestRawWindows numFrames
            |> Map.ofSeq

        let comparisonBestWindows =
            comparisonWindowMap.AllWindowsSequence
            |> WindowMap.ChooseBestRawWindows numFrames
            |> Map.ofSeq

        // merge the two data sets
        let result =
            referenceBestWindows
            |> Map.toSeq
            |> Seq.map (fun (nindex, referenceData) ->
                    if comparisonBestWindows.ContainsKey nindex then
                        Some(nindex, referenceData, comparisonBestWindows.[nindex])
                    else
                        None
                )
            |> Seq.choose id
            |> Seq.map (fun (nindex, (offsets1, spikeCounts1, _), (offsets2, spikeCounts2, _)) ->
                    seq {
                        // NO meta, then WITH meta
                        yield sprintf "%d %d\t%s\t%s" nindex 0 (System.String.Join(" ", offsets1)) (System.String.Join(" ", spikeCounts1))
                        yield sprintf "%d %d\t%s\t%s" nindex 1 (System.String.Join(" ", offsets2)) (System.String.Join(" ", spikeCounts2))
                    }
                )
            |> Seq.concat
            |> Seq.toList

        let header = "nindex withMeta offsets spikeCounts"
        File.WriteAllLines(filePath, header :: result)

    /// Save the profile comparison data generated by CompareProfileWindows as a comparison of counts in the central mass of the peak data
    static member SaveCentralMassData(filePath, windowSize, numFrames, referenceWindowMap:WindowMap, comparisonWindowMap:WindowMap) =

        let referenceBestWindows =
            referenceWindowMap.AllWindowsSequence
            |> WindowMap.ChooseBestRawWindows numFrames
            |> Map.ofSeq

        let comparisonMap = Map.ofSeq comparisonWindowMap.AllWindowsSequence

        // merge the two data sets
        let result =
            referenceBestWindows
            |> Map.toSeq
            |> Seq.map (fun (nindex, referenceData) ->
                    if comparisonMap.ContainsKey nindex then
                        let overlaps =
                            comparisonMap.[nindex]
                            |> WindowMap.SelectOverlappingWindows referenceData

                        if overlaps.Length > 0 then
                            let comparisonData =
                                if overlaps.Length > 1 then
                                    overlaps
                                    |> WindowMap.ChooseBestRawWindow numFrames
                                else
                                    overlaps.Head
                            Some(nindex, referenceData, comparisonData)
                        else
                            printfn "WARNING: no overlapping windows found in comparison set for %d" nindex
                            None
                    else
                        None
                )
            |> Seq.choose id
            |> Seq.map (fun (nindex, (_, spikeCounts1, totalSpikes1), (_, spikeCounts2, totalSpikes2)) ->
                    seq {
                        let centralMassCount1 = WindowMap.CountCentralMass windowSize spikeCounts1
                        let centralMassCount2 = WindowMap.CountCentralMass windowSize spikeCounts2

                        // NO meta, then WITH meta
                        yield sprintf "%d %d %d %f" nindex 0 centralMassCount1 (float centralMassCount1 / float totalSpikes1)
                        yield sprintf "%d %d %d %f" nindex 1 centralMassCount2 (float centralMassCount2 / float totalSpikes2)
                    }
                )
            |> Seq.concat
            |> Seq.toList

        let header = "nindex withMeta centralMassCount ratio"
        File.WriteAllLines(filePath, header :: result)
